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ABSTRACT 

We consider the spectral interpolation error for both Chebyshev pseudo- 
spectral and Galerkin approximations. We develop a family of functionals 
I (u), with the property that the maximum norm of the error is bounded by 
I r (u)/J r , where r is an integer and J is the degree of the polynomial 
approximation. These functionals are used in an adaptive procedure whereby 
the problem is dynamically transformed to minimize I r (u). The number of 
collocation points is then chosen to maintain a prescribed error bound. The 
method is illustrated by various examples from combustion problems in one and 
two dimensions. 
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I. INTRODUCTION 


The objective of this paper is to describe an adaptive pseudo-spectral 
method for solving systems of reaction, diffusion, convection equations. We 
introduce a family of functionals which appear to accurately estimate the 
spectral interpolation error. These functionals are used in an adaptive pro- 
cedure in which the equation is transformed into a coordinate system In which 
the functional is minimized. The number of collocation points is dynamically 
modified to maintain a constant error bound. 

In [1], an adaptive procedure was introduced in which a coordinate trans- 
formation was chosen to minimize the weighted second Sobolev norm of the solu- 
tion. A similar adaptive procedure was presented in [6]. Our numerical 
computations indicate that the new functionals introduced here provide a con- 
siderably sharper error estimate and thereby permit a more efficient solution 
of the differential equations. 

The numerical procedure is illustrated with examples of unsteady axisym- 
metric as well as steady non-axisymmetric flames occurring in gaseous combus- 
tion. This problem is described by the diffusional thermal model in which the 
thermal expansion of the gas is neglected [9], We further simplify by assum- 
ing that the reaction is controlled by a single deficient component so that 
only one species needs to be accounted for, and assume global one step 
Arrhenius kinetics. Thus the model consists of a system of reaction, diffu- 
sion, convection equations for the temperature and the concentration of the 
deficient component which limits the reaction. 

Because of the reaction term, the system has an exponential nonlineari- 
ty. Spurious numerical oscillations can seriously degrade the accuracy of the 
computed solution, primarily through the exponentially nonlinear reaction 



- 2 - 


term. Generally, the reaction term is important only in a narrow region, the 
reaction zone, in which the solution exhibits rapid changes. The Chebyshev 
pseudo-spectral method is subject to oscillations in regions of large gradi- 
ents [5]. We have found that the adaptive procedure described here can be 
effective in reducing these oscillations. An adaptive procedure for a finite 
difference solution of problems of this type is given in [11]. 

This paper is organized as follows. In Section 2, we describe a specific 
model and illustrate the features of the solution which make its numerical 
computation difficult. In Section 3, we describe the adaptive procedure, and 
in Section 4 we describe the new functionals which are used to estimate the 
error. In Section 5, we illustrate the effectiveness of the procedure with 
examples from the model problem described in Section 2. 


II. PROBLEM DESCRIPTION 

To illustrate the behavior of the solution in typical combustion prob- 
lems, we study the following problem for a flame stabilized by a line source 

of fuel. We consider the nondimensional reaction, diffusion, convection sys- 

T _T u 

tem for the reduced temperature 0 = j and reactant concentration C, 

b u 

where T u and T^ represent the temperature in the unburned fresh mixture at 
the source and the burned temperature at infinity, respectively. In polar co- 
ordinates r and <j> , the nondimensional equations are 


© t = A9 — + CN exp(N(l-a)(0-l)/(o+(l-cr)0) 

kC 2 2 

C t = ££ - - CN exp(N(l-a)(0-l)/(a+(l-a )0). 


( 2 . 1 ) 



-3- 


Here o = T^/T^, A is the Laplacian, N is the non-dimensional activation 
energy, and the terms on the right side of (1), respectively represent diffu- 
sion, convection, and a global one-step Arrhenius reaction. The boundary con- 
ditions are given by 

0 = 0, C = 1 for r = 0; 0+1, C+0 for r + °°. (2.2) 


The parameters in (2.1) are the Lewis number, L, the strength k of the 
source of fuel, N, and a. To limit the size of the computational domain, we 
impose the boundary conditions at fixed points rj, r2» rather than at 0 
and 00 . We have verified that the solutions of the problems presented here 
are insensitive to the positions of the artificial boundaries. 

The system (2.1)-(2.2) was analyzed in [8]. It was shown that in the 
limit M = N(1 - a) + «° there exists an axisymmetric solution of the form 


0 - (£>* + 0 ( 1 ) 

C = (1 - 0) + o(i). 


(2.3) 


Furthermore, there exists a critical value of L, L < 1, such that for 

C 1 

L < L the solution (2.3) is unstable to angular perturbations, and arbi- 
C i 

trary initial data close to (2.3), evolve into stationary cells. Computations 
of such cells are presented in [2,3]. 

Another regime of diffusional thermal Instabilities occurs for values 
of L > 1. Based on the analysis by M. Garbey, et al.*, there exists another 


X M. Garbey, G. K. Leaf, H. G. Kaper, and B. J. Matkowsky. in preps*' -it i or. 
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critical value of L, L > 1, such that for L > L the solution (2.3) 

c 2 c 2 

loses its stability to a time oscillatory solution via a Hopf bifurcation. We 

have found this axisymmetric Hopf bifurcation for (2.3) numerically. 

Increasing L along this bifurcation branch, we find that the sinusoidal 

pulsations which occur for L near L become progressively sharper and 

C 2 

steeper and take on the character of relaxation oscillations. At another 

critical value L , a period doubling secondary bifurcation is identified. 

C 3 

The adaptive procedure developed here will be illustrated by a typical 
case after the period doubling secondary bifurcation. Specifically, we con- 
sider the case L = 13, a = .5, N = 40, and k = 11.7. More physically 

realistic values of L occur for other values of the parameter a, N, and 
k. In Figure 1, we illustrate the temporal behavior of 0 and C at the 

fixed point r = 27.4 over a time interval slightly greater than one peri- 

od. The sharp spikes of alternating height are apparent. Computing the solu- 
tion during the spikes consumes most of the cost of the computation. 

To see how the spikes relate to the spatial structure of the solution, we 

exhibit in Figures 2a and 2b the spatial profiles of 0 and C at the four 
times indicated in Figure 1. It is apparent that for all of the time points, 
the solution exhibits very rapid changes over a small spatial interval. Note 
that the concentration profile has a more rapid spatial variation than the 
temperature profile, since L > 1. It is apparent from Figure 2a and 2b that 
the solution becomes much sharper at the time when the temporal spikes 
occur. In the limit M -»■ «°, both 0 and C have discontinuous first 
derivatives. The point of discontinuity is called the flame front, and the 
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reaction term becomes a surface delta function on the front, the strength of 
which was derived in [9], The analysis in [8] is based on this limit. For 
M large but finite, the reaction term is important only in a thin region 
called the reaction zone. 

In computing the solution with the Chebyshev pseudo-spectral method, in- 
accurate resolution of the region of rapid variation can lead to spurious 
oscillations [5,12]. These oscillations can affect the reaction term and give 
rise to a variety of errors, including computed solutions which appear to be 
nonperiodic or chaotic. If the oscillations become large enough, negative 
concentrations occur, which lead to numerical instability. 

To efficiently compute solutions using the Chebyshev pseudo-spectral 
method, we adaptively transform the coordinate system so that in the new co- 
ordinate system the solution appears to vary more gradually. In this way the 
solution can be approximated with a smaller number of collocation points. In 
addition, it is desirable to dynamically add and subtract collocation points, 
as the resolution requirements are more severe when the solution spikes. Our 
procedure is to choose the number of collocation points to maintain a. constant 
error bound. 

These constraints give rise to the need for a good definition and esti- 
mate of the errors associated with the pseudo-spectral method. In the next 
two sections, we describe the adaptive method in general and then discuss the 
functionals we have developed to monitor the spectral interpolation error. 
Although this analysis is only for the error in approximating a function as a 
sura of Chebyshev polynomials, the numerical examples indicate that the pro- 
cedure appears to control the total pseudo-spectral error in solving the equa- 


tions of the mathematical model. 
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III. ADAPTIVE PROCEDURE 

We consider, for purposes of discussion, the model reaction, diffusion, 
convection equation in one space dimension 

u fc = u xx + u x + R(u), -i < x < 1, (3.1) 

where R(u), is a nonlinear function. The initial and boundary conditions are 
not important for the purpose of this discussion. In the Chebyshev pseudo- 
spectral method, the solution to (3.1) is approximated as a finite sum of 
Chebyshev polynomials 

J 

u = l a (t)T (x) , (3.2) 

j=0 J J 

where Tj(x) = cos(jcos ^x) and the coefficients aj are obtained from 

collocating the solution at the points 

{Xj = cos(^-), j = 0, •••,!}. (3.3) 

For a discussion of this method, see [5]. The method differs from finite 
differences in that it is spectrally accurate; i.e., for sufficiently smooth 
solutions the error decays faster than J -r for any power of r. In partic- 
ular for analytic functions, the error decays exponentially [13]. A precise 
formulation of this property is discussed below. 

It Is known that this method can exhibit oscillations when solutions with 
rapid spatial changes, such as those exhibited in Figures 2a, b, are 

computed. In order to improve the ability of the method to accurately compute 
such solutions, an adaptive procedure was introduced in [1]. 
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In this procedure, we dynamically vary the coordinate system for (3.1). 
We transform (3.1) by the mapping 


x = q(s,a), 


where s is the new independent variable. Here q is a given function, 
and a is an unknown parameter vector to be determined. In [1] and in the 
computations presented below, the function q was taken to be 


q(s,a) = — tan *[a. tan 7- (s' - 1)] + 1, 

~ IT 14 


where 


a„-s 

s' = r 

a 2 s_1 


a. > 0, -1 < ct 9 < 1 


and a = (a^^). 

To determine a, we suppose that I r (u) Is a functional of the solu- 
tion and that some measure of the error is bounded by I r (u)/jr for some 
integer r. Then I f (u(q(» ,a) ) ) Is a measure of the error in the coordi- 
nate system determined by a . We can then minimize the upper bound on the 

error by choosing a such that I r (u(q,a)) = T r (a) is minimized. If e 
is a desired error level, then J can be determined as the smallest integer 
such that 


1 (a) 

j r > J L ~ . 


(3.4) 


In this procedure, it is necessary to interpolate the solution to the new 
collocation points x^ = q(sj,a), where (sj) are the Chebyshev colloca- 
tion points (3.3). This can be accomplished without loss of spectral accuracy 
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by evaluating the global approximation (3.2). From this discussion it is 
clear that the role of l r (u) is crucial. In [1,6], I r (u) was chosen as 
the weighted second Sobolev norm 

I 2 (u) = [/ dxw(x) [ | u xx | 2 + |uj 2 + |u| 2 ]] 1/2 , 
where w(x) is the Chebyshev weight function 

w(x) = (1 - x 2 )- 1/2 . 


The adaptive procedure is therefore to find a such that 


1 


1/2 


Ka) -J ds[|u gs | + [uj + |u| ]w(s) 


(3.5) 


is minimized. 

The choice of (3.5) is motivated by the estimates obtained in [4], where 
it is shown that if u is smooth, and PjU the Chebyshev interpolant, then 
the error ej = u - Pju satisfies 


lull. 


le j"o < C r 


(3.7) 


for J > 0. Here II II q is the weighted L 2 norm, 

II e jll q = / ds w(s) hjl 2 , 
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and II II is the weighted Sobolev norm. We make the following re- 
marks about (3.7) 

a) As is typical in pseudo-spectral methods, (3.7) is a family of error 
bounds and the error seen in practice is the smallest of the right hand sides 
of (3.7). Thus, any adaptive procedure can only minimize an upper bound on 
the error. For the combustion problems discussed here, we choose r = 2, i.e. 
(3.5), although when the spatial profiles become more smoothly varying, this 
may overestimate J. 

b) The error bound on the left hand side of (3.7), which is to be bounded, 
depends on the particular coordinate system that is chosen and may not neces- 
sarily relate to the quality of the computed solution. 

c) The estimates leading to (3.7) are not sharp. 

d) It has been observed that (3.7) does not yield good coordinate systems 
for r = 1 [1,6] . 

e) In practice we find that (3.5) is not sufficiently sensitive to variations 
in the behavior of the solution, when compared with the functional described 
below. The use of (3.5) tends to overestimate the number of collocation 
points J. A discussion of a new family of functionals, which appear to 
better track the spectral interpolation error, now follows. 
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IV. THE NEW FUNCTIONAL 

To better understand the behavior of the error in the Chebyshev method, 
and to get a better error functional, we start by considering a smooth func- 
tion 

1 f(x)T (x) 

* dx, 

1 

where cq = 2; c^ = 1 for k > 0, and its Chebyshev-Galerkin approximation 

f j-i (x) ' p j-i e(x) ■ J„ a kV x> - 

k=0 

We are interested in an error estimate that will not depend on the coordinate 
system; otherwise it will be impossible to know whether the new coordinate 
s is really better than the original one. The most natural norm in which to 
measure the error is the maximum norm. Thus we look at the pointwise error 


f(x) = l a,T (x); 
k=0 K K 


- 7 -/ 


“k - 


|f(x) - fj_j(x)| = | l a k T k (x)| 

k=J 


(4.1) 


Consider first the leading term, aj» in (4.1). In fact, since for k > 0, 




x V' - - 


we obtain 


k 2 T, 


fa 


= ~ — — ^ J 1 (A-x 2 T')' f (x) dx. 

Cj J -1 


(4.2) 


We then integrate by parts, noting that the boundary terms vanish, to obtain 
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a J = 


V 


/ /l— X 2 T' f'(x)dx = - / (/l-x 2 dx, (4.3) 


C J J - 1 


or upon defining L = /l-x 2 , 


= - -^r / a 2 f) 


dx. 


V 


/l-x 2 


(4.4) 


This formula is a special case of the formula derived in [5] for a general 
Sturm-Liouville problem. We repeat the integration by parts to obtain 


a J = f (L2kf) ~~~ dx » k " i* 2 , 


(4.5) 


C J J 


AZ 2 


From (4.4), we can bound aj by using the Cauchy-Schwartz inequality 


a 2 f ) 2 dx 1 i/2. 1 rr Tj 
' J '” j7LJ AA?" 


•I <ly[/ 


At? 


dx] 


1/2 


(4.6) 


The last factor on the right hand side of (4.6) is one, by the normalization 
of Tj . It is therefore reasonable to take 


1 ,. 2_2 i 


! [f] = {/ (L f) dx} 2 = {/ J*— u/iA 2 ^) 2 f] 2 f 


~ l /l-x 2 




(4.7) 


instead of (3.5). Similarly, we can define 


i 2k [ f , .[/ 


_1 /l-x 2 


and take the minimum over k as a measure for the error. Going back to 
(4.1), we can bound the (Galerkin) spectral interpolation error in the maximum 
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norm 


I,(f) I A (f) 

max | f(x) - f . , (x) | _< l | a, | < — 5 — + — * — • (4. 8) 

-i<x<i k=j K r 


Several comments are now in order: 


a) The functional (4.7) gives different weight to the boundaries than (3.5) 
and weights the second derivative differently from the first derivative. 


b) The same analysis, though integrating by parts only once, shows that 




dx, 


where K is some constant. Thus even the lowest order error bounded does not 
reduce to the weighted Sobolev norm. This may explain the lack of success of 
using the weighted Sobolev norm of the first derivative in [1,6]. 


c) The bound (4.7) can be extended to bound the spectral interpolation error 
based on collocation, with the loss of at most a factor of log J. This is 
because the right hand side of (4.8) is also an upper bound for the best poly- 
nomial approximation to f in the sup norm. It is known [10] that the max- 
imum norm of the error of the interpolant of f, at the points (3.3) will dif- 
fer from the error obtained from the minimax polynomial, by a factor that 
grows at most logarithmically. 
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d) The analysis presented here is only for the spectral interpolation 
error. However, the numerical results indicate that minimizing (4.7) is 
effective in controlling the total error in our computations. 


V. NUMERICAL RESULTS 

In this section we present some numerical results illustrating the ef- 
fectiveness of (4.7) in tracking the spectral interpolation error. We con- 
sider the model described by (2.1) and the parameters described in Section 
2. In our adaptive procedure we consider the current coordinate system 


described by 


a = a 

/V A/ 


The solution is updated in time and I 2 (u) is 


computed at every time step, 
triggered whenever 


A search for a new coordinate system is 


or 


I 2 ( u ) 

*min 


I 2 (u) 

^min 


< Cp 


y q~2 * 


(5.1a) 


(5.1b) 


where I ^ is the value of l2(u) at the previous time that a coordinate 
system was found. We observe that (5.1a) may indicate that the number of col- 
location points should be reduced in order to more efficiently compute the 
solution, whereas (5.1b) may indicate that more collocation points are re- 
quired to achieve the desired accuracy. Both (5.1a) and (5.1b) indicate that 
the solution has changed sufficiently from the previous time a coordinate sys- 
tem was found, so that a new minimization should be attempted. For the compu- 
tations used below, we took cj “ .3, c 2 ■ 1.7. At present only one coordi- 
nate system is used for both variables. The function actually used in the 
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evaluation of ^(u) is a linear combination of 0 and C. For the compu- 
tations of pulsating flames presented below, we use u = .30 + .7C. In ad- 

dition, we specify a maximum number of collocation points, J max , which is to 
be used in case the program cannot meet the specified error. 

In Figures 3a and 3b, we display the temporal behavior of 0 and C 
for the same computation as shown in Figure 1. In these figures, we also plot 
the number of collocation points as t varies. In Figure 3a, the choice of 
a is obtained from (4.7), while in Figure 3b the choice of a was ob- 
tained from (3.5). The value of e in each case was 0.0025. 

It is apparent that the use of (4.7) is more sensitive to the behavior of 
the solution than the use of (3.5), and permits a more rapid reduction in the 
number of collocation points near the spike, where the maximum cost of the 
computation occurs. In fact, the computation for Figure 3a required roughly 
30% of the time of the computation for Figure 3b. 

In Figures 4a and 4b, we plot the spatial profiles for 0 and C at 
the initial spike (t ~ 2.9). In each case we took the spatial profile at 
the time that the functional was largest. The value of the square of the 
functional is indicted in the figures. The major source of spatial errors is 
the oscillations in the profiles, primarily in the concentration profile. We 
observe that the solution in Figure 4b (using (3.5)) is slightly more oscilla- 
tory than the solution in Figure 4a (using (4.7)), although both solutions are 
acceptable. The effect of increasing the error tolerance is shown in Figures 
5a (using (4.7)) and 5b (using (3.5)). In this case we set e = .0035. It 
can be seen that Figure 5a is more accurate. In this case the maximum number 
of collocation points for Figure 5a was 321 while in Figure 5b the error level 
was not met with 361 points (the maximum number of points allowed in this run 
of the program). The ratio of the running time was 5 to 1. 
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As another example, we reduced the maximum allowed number of collocation 
points to 241. Spatial profiles using (4.7) (Figure 6a) and (3.5) (Figure 6b) 
are illustrated. Although both solutions are degraded, it is apparent that 
Figure 6a is considerably more accurate than Figure 6b. 

We next discuss the application of the method described above to the 
computation of stationary cellular flames. We consider (2.1), allowing both 
radial and angular perturbations, for values of L < 1. We consider L < L 

C 1 

so that the axisymmetric solution to (2.1) is unstable to angular perturba- 
tions and arbitrary initial data sufficiently close to the axisymmetric solu- 
tion, will evolve to a stable stationary cellular solution. These cellular 
flames, sometimes called wrinkled flames, are characterized by the appearance 
of pointed crests and troughs in the reaction zone. The crests point in the 
direction of the burned products of combustion and are cooler than the 
troughs. This behavior, described analytically by the interaction between the 
fundamental and its first harmonic [8], is also observed experimentally [7] 
and is present in our computations at finite activation energy as well [2,3]. 

In our computations, we employ a Fourier pseudo-spectral method in the 
angular direction <}> and a Chebyshev pseudo-spectral method in the radial 
direction r. At present we employ a transformation in the radial direction 
with a independent of <J> , and therefore we require sufficient radial 
resolution to resolve the reaction zone in each angular direction. The sta- 


tionary state is achieved by integrating in time until the solution 

3 0 3 C 

equilibrates. The convergence is monitored by the maximum of |g-^| , | _ g^| » 
over all collocation points. We require this maximum to be less than 10 - ^. 


The temporal integration uses the backward Euler scheme with operator split- 


ting. The reaction terms are treated explicitly while all other terms are 


treated implicitly [2], 
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We illustrate our method with a computation of a three cell. The 
parameter values employed are k = 11.0, L = 0.44, N = 20, and a - 0.615. 
We first describe the behavior of the solution. Figures 7a, b illustrate the 
temperature and concentration, as a function of (j> , at two different radial 
positions. We observe that for r = 9.34, the temperature exhibits a large, 
nearly sinusoidal oscillation with few visible higher frequencies. For r = 
10.15 the angular variation is reduced, but the temperature exhibits crests 
and troughs. This is due to the interaction of the fundamental and its first 
harmonic; the other harmonics are very small. This effect is illustrated in 
Figure 7c by plotting the mean, fundamental, and first harmonic as a function 
of r. 

The data in Figure 7 were obtained by using 101 points in r and 128 
collocation points in <j> . This solution has been validated by further grid 
refinements, and for this number of radial collocation points, the solution 
does not depend upon the specific functional used to determine the coordinate 
system. To illustrate the effectiveness of (4.7), we consider a cruder 
approximation, using a reduced number of collocation points. 

We compare the fine grid solution with computations using (4.7) and (3.5) 
with 49 collocation points. The objective function is obtained by averaging 
the functional over all angles. In Figure 8a we plot the mean value of the 
solution, comparing the fine grid solution with the coarse grid solution ob- 
tained by using (4.7). In Figure 8b we do the same for the fundamental and 
its first harmonic. In Figures 9a and 9b we present the same data, except 
that the coarse grid solution is obtained by using (3.5). 

It is apparent from the figures that the coarse grid solutions are con- 
siderably less accurate on the angular variations than on the mean. Although 
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the first harmonic is small relative to the mean and the fundamental, we con- 
sider it an important part of the computation, as the pointed crests and 
troughs in the reaction zone are features predicted by analysis and by our 
computations and are observed experimentally. It is apparent from the figures 
that by using (4.7), we obtain a more accurate solution than by using (3.5). 
A more detailed numerical study of cellular flames is presented in [2,3]. 


VI. CONCLUSION 

We have presented a functional that appears to provide a relatively sharp 
upper bound for the spectral interpolation error, in the maximum norm, for the 
types of solutions occurring in typical combustion problems. The functional 
is derived directly from the Chebyshev expansion of a smooth function; it is 
the first member of a family of functionals, depending upon successively high- 
er derivatives, which bound the spectral interpolation error. We have em- 
ployed the functional to both reduce the error for a fixed number of colloca- 
tion points and to dynamically vary the number of collocation points, depend- 
ing on the magnitude of the functional. The adaptive procedure based on this 
functional is shown to be more effective than the procedure based on the 
weighted Sobolev norm, in computations of both pulsating and cellular flames. 
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FIGURE CAPTIONS 

1. 0 and C as a function of time at a fixed radial location. 

2a. Spatial profiles of 0 at four different times. 

2b. Spatial profiles of C at four different times. 

3a. 0, C, and J as a function of time. J and a were obtained from 

minimizing (4.7). 

3b. 0, C, and J as a function of time. J and a were obtained from 

rsj 

minimizing (3.5) 

4a. Spatial profiles of 0 and C at the time when the objective func- 
tional is largest. J and a were obtained from minimizing (4.7) 

and e = .0025. 

4b. Spatial profiles of 0 and C at the time when the objective func- 
tional is largest. J and a were obtained from minimizing (3.5) 

and e = .0025. 

5a. Spatial profiles of 0 and C at the time when the objective func- 
tional is largest. J and a were obtained from minimizing (4.7) 


and e 


.0035 
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5b. Spatial profiles of 0 and C at the time when the objective func- 
tional is largest. J and a were obtained from minimizing (3.5) 

and e = .0035. 

6a. Spatial profiles of 0 and C at the time when the objective func- 
tional is largest. J and a were obtained from minimizing (4.7), 

with J constrained to be _< 241. 

6b. Spatial profiles of 0 and C at the time when the objective func- 
tional is largest. J and a were obtained from minimizing (3.5), 

with J constrained to be 241. 

7a. Angular variation of 0 and C for a cellular flame; r = 9.34. 

7b. Angular variation of 0 and C for a cellular flame; r = 10.15. 

7c. Mean, fundamental and harmonic for a cellular flame, as a function of r. 

8a. Mean of 0 as a function of r. Coarse grid solution obtained from 

minimizing (4.7). 

8b. Fundamental and harmonic of 0 as a function of r. Coarse grid solu- 
tion obtained from minimizing (4.7). 

9a. Mean of 0 as a function of r. Coarse grid solution obtained from 
minimizing (3.5). 
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9b. Fundamental and harmonic of 0 as a function of r. Coarse grid solu- 
tion obtained from minimizing (3.5). 
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Figure 2b. 
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Figure 5a. 
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Figure 5b . 
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Figure 9a. 
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